Collisionally inhomogeneous Bose-Einstein condensates in double-well potentials 
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In this work, we consider quasi-one-dimensional Bose-Einstein condensates (BECs), with spatially 
varying colhsional interactions, trapped in double well potentials. In particular, we study a setup in 
which such a "collisionally inhomogeneous" BEC has the same (attractive-attractive or repulsive- 
repulsive) or different (attractive-repulsive) type of interparticle interactions. Our analysis is based 
on the continuation of the symmetric ground state and anti-symmetric first excited state of the non- 
interacting (linear) limit into their nonlinear counterparts. The colhsional inhomogeneity produces 
a saddle-node bifurcation scenario between two additional solution branches; as the inhomogeneity 
' becomes stronger, the turning point of the saddle-node tends to infinity and eventually only the 

two original branches remain present, which is completely different from the standard double- well 
^ , phenomenology. Finally, one of these branches changes its monotonicity as a function of the chemical 

O ■ potential, a feature especially prominent, when the sign of the nonlinearity changes between the two 

' wells. Our theoretical predictions, are in excellent agreement with the numerical results. 

a^ ; 

I ^' I. INTRODUCTION 

c/^ : 

' The remarkable progress in the experimental and theoretical studies of Bose-Einstein condensates (BECs) P, Q 
■ over the past two decades has sparked an intense study of the coherent nonlinear structures that arise in this setting. 
' Such structures are rather naturally expected to emerge due to a well-established mean-field description of BECs 
^ by the Gross-Pitaevskii (CP) equation, namely, a nonlinear Schrodinger (NLS) equation, in which the nonlinearity 
is introduced by the interatomic interactions. This feature has inspired many relevant theoretical and experimental 
. studies devoted to macroscopic nonlinear matter waves, such as bright matter- wave solitons in attractive BE Cs 
J> ' as well as dark 0, H, and gap [TO] matter wave solitons in repulsive BECs (see also the recent review [Tlj|). One 
Tj" , of the important elements of the versatility of this atomic physics setting is the existence of diverse types of external 
' potentials in the CP equation accounting for the electric, magnetic, optical or combined confinement of ultracold 
' dilute alkali vapors that constitute the BEC. Some of the most typical forms of such trapping potentials include a 
. I parabolic and a spatially periodic one (created by the interference of counter-propagating laser beams, the so-called 
optical lattice). The NLS equations with similar potentials are also relevant in the context of nonlinear optics, where 
they model the evolution of optical beams in graded-index waveguides and periodic waveguiding arrays [12l [l3| . 

One interesting type of potential that has drawn a considerable amount of attention is the double-well potential. 
This may originate, for instance, from the combination of a parabolic with a periodic potential. Its experimental 
^ ■ realization was featured in Ref. [14], where a variety of interesting phenomena were studied; these include Josephson 
oscillations and tunneling for a small number of atoms, or macroscopic quantum self-trapping and an asymmetric 
partition of the atoms between the wells for sufficiently large numbers of atoms. In parallel to these experimental 
findings, numerous theoretical insights on this topic have emerged [15, 16, 17, jl8., JJ,, 20, 21, 22, 23]. These concerned 
finite-mode reductions, analytical results for specially designed shapes of the potential, quantum depletion effects, 
and other theoretical aspects. Interestingly, double-well potentials have also been studied in app lications arising in 
the context of nonlinear optics, including twin-core self-guided laser beams in Kerr media [24| . optically induced 
dual-core waveguiding structures in a photorefractive crystal (25| , trapped light beams in a structured annular core 
of an optical fiber [26[, and so on. It is relevant to point out here that double well settings have been examined 
not only in one-component systems, but also in multi-component cases. In particular, a recent study motivated by 
two-component BECs can be found in [13], while similar attempts have been made in the context of the so-called 
spinor BECs (where it is possible to have three, and even five components) in the works of [H, [2^. These works 
examined not only finite-mode reductions of the multi-component case, but also phenomena beyond the level of the 
mean-field description such as quantum entanglement and spin-squeezing properties. 

On the other hand, nonlinear matter-waves have not only been studied in a variety of external potentials, but also 
in the presence of temporally or spatially varying external fields manipulating the interatomic interactions. Indeed, 
the s-wave scattering length (which characterizes the nonlinearity coefficient in the CP equation) can be adjusted 



00 

o 



- '—I 

X 



experimentally by employing either magnetic [3Ctl3lj or optical Feshbach resonances [32| in a very broad range. This 
fiexibility of manipulation of the interatomic interactions has motivated a significant number of studies both on the 
theoretical and on the experimental front; in particular, on the experimental side, the formation of bright matter-wave 
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solitons and soliton trains for ^Li [1,!3] and ®^Rb ^ atoms used a tuning of the interatomic interactions from repulsive 
to attractive. Also, this type of manipulations was instrumental in achieving the formation of molecular condensates 
[3^ . and the probing of the BEC-BCS crossover [s^l- On the other hand, theoretical studies have predicted that 
a time-dependent modulation of the scattering length can be used to stabilize attractive higher-dimensional BECs 
against collapse ^35], or to create robust matter- wave breathers in lower-dimensional BECs [36J- While the above 
studies focused on temporal variations of the interaction strength, more recently spatial variations of the nonlinearity 
have come to be of interest in the so-called "collisionally inhomogeneous" environments. These have been found to 
lead to a variety of interesting developments including (but not limited to) adiabatic compression of matter-waves 
[13, [13 1 Bloch oscillations of matter- wave solitons [33, atomic soliton emission and atom lasers [1^, enhancement 
of transmittivity of matter- waves through barriers [io, [4l| , dynamical trappin g of matter- wave solitons [i^l , stable 
condensates exhibiting both attractive and repulsive interatomic interactions [42l |. the delocalization transition of 
matter waves )43|. and so on. Many different types of spatial variations of the nonlinearity have been considered, 
including linear |37l. lioj . parabolic [ill, random [4^1 , periodic [isl. l46l . l47j . and locahzed (step-like) [H, [H, [4§| ones. 
Furthermore, a number of detailed mathematical studies [13, ISlL l52l| have appeared, addressing aspects such as the 
effect of a "nonlinear lattice potential" (i.e., a spatially periodic nonlinearity) on the stability of matter- wave solitons, 
and the interplay between drift and diffraction/blow-up instabilities. More recently, the interplay of nonlinear and 
linear potentials has been examined in both continuum [s^ and discrete settings. 

The aim of the present work is, in fact, to combine these two interesting settings, namely the double-well potential 
and a collisionally inhomogeneous (i.e., spatially dependent) nonlinearity. As, arguably, one of the simplest forms of 
this combination, we select a coefficient of the nonlinearity which is piecewise constant, in line with the suggestions 
of [3^ |49| . [55! and can vary between two (smoothly connected) pieces with the same sign, or even two pieces with 
opposite signs. The latter form of a spatially dependent nonlinearity appears to be one that should be straightforwardly 
experimentally realizable 15611 , through the use of magnetic field gradients of moderate size for atom chips (see also 
the relevant discussion in [55|). The phenomenology that we observe in this setting appears to be remarkably different 
from that of the standard double- well. In particular, as we detailed in an earlier relevant work [23], the linear states 
of the underlying potential can be continued into nonlinear ones in the presence of the pertinent cubic nonlinearity. 
As prototypical examples of these linear states one can consider the symmetric ground state and the anti-symmetric 
first excited state of the double-well potential. A remarkable feature, however, which takes place for sufficiently 
strong nonlinearity (to the symmetric state for attractive interactions and to the anti-symmetric state for repulsive 
interactions) is a symmetry-breaking bifurcation of the pitchfork type. This bifurcation generates two new asymmetric 
solutions (which exist for sufficiently large nonlinearity). What we find here is that even a weak inhomogeneity in 
the collisional properties of the condensate changes the nature of this bifurcation from a pitchfork to a saddle-node 
one. This may be anticipated given the "non-parametrically robust" nature of the pitchfork bifurcation which will 
yield similar changes in the presence of different asymmetries (see also [l^). In addition to that, increasing the 
inhomogeneity strength shifts (eventually to infinity) the turning point of the saddle-node bifurcation. Thus, for 
sufficiently large variation of the nonlinear coefficient between the two wells, only two nonlinear states emanate from 
the ground and first excited states of the linear problem. Furthermore, for one of these states, the monotonicity 
of the number of particles changes as a function of its nonlinear eigenvalue parameter (i.e., the chemical potential). 
This is a rather unusual feature that leads to a bifurcation diagram entirely different from those of the homogeneous 
interatomic interactions (i.e., homogeneous nonlinearity) limit of either the one-component [23 | or of the multi- 
component system (see e.g., 127] V 63] These traits are captured accurately, as we will show below, by a two-mode, 
Galerkin-type approximation [2l|, |2^ applied to the present collisionally inhomogeneous double well setting. 

Our presentation is structured as follows. The model and the semi-analytical predictions regarding its stationary 
solutions are presented in section II. Our corresponding directly numerical findings and their comparison with the 
analysis are given in section III. In section IV, we present an alternative theoretical viewpoint more akin to the model 
dynamics and compare its results to the full results of the partial differential equation. Finally, in section V, we 
present our summary and conclusions. 



II. THE MODEL AND ITS SEMI-ANALYTICAL CONSIDERATION 



The dynamics of a quasi one-dimensional (ID) condensate, oriented along the a;-axis, can be described by the 
following GP equation (see, e.g., Ch. I in [ll[ and [13]), 

^h^t^^(^~^ + Vix)+g\^\^-^?j^, (1) 

where ^'(a;, t) is the mean-field order parameter, m is the atomic mass, and fi is the chemical potential of the effectively 
ID system. The nonlinear coefficient g arises from the interatomic interactions and has an effective ID form, namely 
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g = 2hu!±a, where lu± is the transverse confining frequency and a is the three-dimensional (3D) s-wave scattering 
length (the cases a > or a < correspond, respectively, to repulsive or attractive interatomic interactions). The 
external potential V{x) in the GP Eq. ([T]) consists of a regular harmonic trap and a repulsive potential localized at 
the harmonic trap center, namely, 

V(x) = Imujlx^ + ^osech^ ( — | , (2) 

2 \WB J 

where tOx is the longitudinal confining frequency, while Vq and wb are, respectively, the strength and width of the 
localized potential; the latter is in fact a barrier potential that may be created by a blue-detuned laser beam, repelling 
the atoms in the condensate. It is clear that the combination of the harmonic trap and the barrier potential is in fact 
a double-well potential. 

Moreover, we assume that the collisional properties of the condensate are spatially inhomogeneous, i.e., a = a{x), 
with the function a{x) taking different, but smoothly connected, values in the two wells of the external potential. In 
particular, we consider that an external magnetic or optical field (see below) modifies the scattering length of the 
condensate as follows, 

a{x) ^ ao + aitanh (^-^^ , (3) 

where ag and ai are constant values of the condensate's scattering length in the absence and in the presence of the 
external field, respectively, which are smoothly connected through the tanh function (here, W is the spatial length 
scale on which the transition between the two values oq and ai takes place). Apparently, far away from the harmonic 
trap center at a; = (or, in other words, far away from the barrier), the scattering length takes the values a = ao±ai, 
for X —^ ±oo respectively. Such an inhomogeneity of the scattering length may be realized in practice upon employing a 
bias homogeneous field and imposing a steep localized gradient on top of it. In such a case, in a quasi ID configuration 
(as is the case under consideration) this will lead to a constant scattering length with value a = ai in the left well, 
which is followed by a localized change of a, finally ending up with a second value a = a2 in the right well. Notice 
that in our model we choose the function tanh to analytically describe the transition between these different constant 
values of a since there are no ideal steps, but rather close approximations to it. Also, we naturally assume that we are 
relatively close to a Feshbach resonance, in order to easily manipulate the scattering length with the external field. 

Next, measuring time in units of the inverse transverse trapping frequency (jjJ^, length in units of the transverse 
oscillator length l± = ^/K/mujZ, energy in units of fiu}±, and introducing the normalized wavefunction u{x,t) — 
(2|a^|)i/2*(x,t) (where ar is a reference value of the scattering length) , we reduce the GP Eq. ([1]) to the following 
dimensionless form: 

idtu:^ Hu + r{x)\u\'^u- fiu. (4) 

In the above equation, 

H^-ldl + V{x) (5) 
is the "single-particle" operator in the normalized external potential 

V{x) = -n^x^ + Vq sech" (-) , (6) 
2 \w/ 

with Q = LUx/<-^i_ being the normalized harmonic trap strength and w = 'Wb/1±- The nonlinearity coefficient T(x) in 
Eq. dH) is given by 

r(a;) =a + /3tanh(6a;), (7) 

where a = ai/\ar\, P = a2/\ar\, and b = lj_/W. Apparently, the cases F > or P < correspond to repulsive or 
attractive interactions, respectively. We finally mention that the number of particles in the condensate M is now 
defined as J\f= {l±/2\ar\)N, where 

/ + 00 
\u{x,t)\^dx (8) 

is an integral of motion (normalized mumber of particles) of the GP Eq. Also, in our analysis (and particularly 
in the typical numerical results that will be presented below) we will assume fixed parameter values — 0.1, Vq — 1, 
w = 0.5 and b = 1. 
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In the non-interacting limit [i.e., for r(x) = 0], the spectrum of the underlying linear Schrodinger equation consists 
of a ground state, Uo{x), and excited states, ui{x) {I > 1). Our analytical investigation of the nonlinear problem (at 
the stationary level) in the weakly nonlinear regime consists of applying a Galerkin-type, two-mode approximation to 
Eq. In particular, we assume that the wavefunction u{x, t) can be decomposed as 

u{x,t) = Cf){t)u(){x) + ci{t)ui{x), (9) 

where co(t) and ci{t) are unknown time-dependent complex prefactors. Substituting Eq. ([9]) into Eq. ([4]) and 
projecting the resulting equation to eigenfunctions uq and Ui, we obtain the following equations: 

ico = {luo- ^^)ca + Ao\cofco + Do{2\cQ\^Cl+clcl)+B{2\clfco + clc*„) + Dl\cl\'^Cl), (10) 
ici = (lui- n)ci+Ai\ci\^ci+Di{2\ci\^co + cl4)+B{2\co\^ci+clcl) + Do\co\^co). (11) 

In these equations, dots denote time derivatives, stars denote complex conjugates, and wo,i are the eigenvalues cor- 
responding to the eigenstates uo,i; in the numerical examples presented herein, wq and wi can be numerically found 
to be 0.133 and 0.156 respectively. We should note in passing here that although wo,i and wo,i are numerically 
obtained through a simple eigensolver of the underlying linear Schrodinger problem, it is also, in principle, possi- 
ble to develop a perturbative approach towards obtaining such eigenvalues and eigenvectors following the work of 
[ssj . rendering our semi-analytical approach developed below more proximal to a fully analytical one. The over- 
lap integrals Aq = J r{x)uQ{x)dx, Ai — J r{x)uf{x)dx, B — ^ T{x)u^{x)u\{x)dx^ Dq — J T{x)uQ{x)ui{x)dx and 
Di — J T{x)uQ{x)uf{x)dx are constants. Notice that uq{x) and ui{x) are both real functions (due to the Hermitian 
nature of the underlying linear Schrodinger problem) and are also orthonormal. In what follows, in the definition of 
r{x) in Eq. (O, we wiU set 

a + f3^-l, a-(3^-l + e, (12) 

with < e < 2, as shown in Fig. [TJ Notice that in the limiting case e = 0, the nonlinearity coefficient takes a constant 
value, r(a;) — —1, and the problem is reduced to the study of a coUisionally homogeneous condensate with attractive 
interatomic interactions in a double well potential; this problem was studied in detail in Ref. ^2^ by means of the 
considered two-mode approximation for both symmetric and asymmetric double-well potentials. As a result of this 
choice and of fixing the rest of the parameters, the above overlap integrals Aq, Ai, B, Dq and Di are functions of e 
and their dependence on this parameter is shown in Fig. [5) 
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FIG. 1: The nonlinearity coefficient T{x) — a + /3ta.nh{bx) with parameters a + f3 — —1, a — /3 = —1 + e, and 6=1. The 
figure displays some examples of the shape of r(a;) for difi'erent values of e (i.e., the inhomogeneity parameter): e = 0.5 (blue 
solid line), e = 1 (green dashed line) and e = 2 (red dashed-dotted line.) 



We now introduce in Eqs. (fTO |) -(fTT |l the amplitude-phase variables cj — pje'^^\ j = 0, 1 (the functions pj and Lpj 
are assumed to be real and are time-dependent) and derive from Eqs. (jlOp and (jlip the equations for po and {p^: 

Po = Bpoplsm{2(p) + {DopIpi + Dipl)sin{ip), (13) 

ifio = ip-LJo)-Aapl-B{2pl+plcos{2ip))~{3DopoPi+Di^)cos{ip), (14) 

Po 

where ip = ipo — ^pi is the relative phase between the two modes. The equations for pi and ipi can directly be obtained 
by interchanging indices 1 and in Eqs. ([TH]) and (fH]) . 
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FIG. 2; The dependence of the overlap integrals Ao, Ai, B, Do and Di on e (for our choice of the rest of the parameters of the 
potential and of r(a;)). 



Focusing on real solutions of Eq. (jH), we consider the steady solutions, i.e., po = (po = [associated to the fixed 
points of the dynamical system of Eqs. (jl3[) - (|14p ]. with if = kn with k an integer. In such a case, the equations for po 
[Eq. (113])] and pi are automatically satisfied, while the equations for ipQ [Eq. p^ ] and ipi are reduced to the following 
algebraic system: 



^3 

{p-u;o)-Aopl-3Bpl-3DoPoPi-Di^ = 0, (15) 

Po 

{p^LUi)-Aipl-?,Bpl-iDipipo-Da— = 0. (16) 

Pi 

The semi-analytical part of our considerations will consist of solving these algebraic conditions for given linear potential 
V{x) and nonlinear spatial dependence r(a::) parameters (as mentioned above) and we will compare the findings from 
this two-mode, Galerkin-type truncation with the full numerical results in what follows. 

It is interesting to briefly address the case where (f> ^ kir (hence n{<j)) 7^ 0). In that case, the stationary form of Eq. 
(fT3)) leads to: 

cos(^)^-^4±^. (17) 

Then it is straightforward to observe that since Dq < Di, for all e, and B < Di for e > 0.995 (from Fig. for such 
values of e, the fraction of the right-hand-side of Eq. (|17|) can be shown to be < —1 and, hence, there are no other 
solutions satisfying Eq. p^ . For e < 0.995 we could not perform a similar semi-analytical calculation. However, we 
confirmed numericalUy the non-existence of additional solutions (to the ones with (j) = kn) for this domain. 



III. NUMERICAL RESULTS 



We begin the exposition of our numerical results by considering < e < 1, in which case T{x) < 0, i.e., attractive 
interatomic interactions. The top panels of Fig. [3] show two examples: e = 0.25 (left) and e — 0.5 (right). Each 
of these panels presents the complete diagram of the numerically generated solutions (blue solid lines correspond to 
stable solutions and red dashed lines to unstable ones) to Eq. ^ and the corresponding analytical predictions (green 
dashed-dotted lines), namely, the results of the two mode approximation obtained by solving Eqs. (fT5)) and (fTB|) . The 
solutions are expressed in terms of the normalized number of atoms N [see Eq. ([8])] as a function of the normalized 
chemical potential /i. The branches are obtained using a fixed-point Newton iteration and numerical continuation 
from the non-interacting limit of the system (i.e., in the absence of nonlincarity). Then, numerical linear stability 
analysis is performed to determine whether each branch does or does not possess real eigenvalues (which result in 
instability). It is readily observed that the analytical predictions for the different branches through Eqs. (|15p - (fT6l) 
are in very good agreement with the numerically found solutions. 

Let us now discuss the various branches appearing in the bifurcation diagrams in more detail. Branch I corresponds 
to the asymmetric (due to the coUisional inhomogeneity) solution starting at p. = luq, the eigenvalue corresponding 
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FIG. 3: Top panels: The normalized number of atoms A'^ [see Eq. (|8])] of the solutions of Eq. Q for the case of attractive 
interatomic interactions i.e., for a nonlinearity coefficient r(a::) with parameters a + (3 — —1 and a ~ l3 — —1 + e, with e = 0.25 
(top left) and e = 0.5 (top right), as a function of the normalized chemical potential /i. The blue solid lines and red dashed 
lines denote the stable and unstable numerically found solutions. The green dashed-dotted lines depict the result of the two- 
mode approximation. Notice that as — > 0, the spatial profiles of the two branches tend to the linear eigenmodes uo,i and 
accordingly jj, — > a;o,i. Bottom panels: The profiles of the wave functions corresponding to branch I (upper blue) and II (lower 
green). Along each branch, the profiles are shown for two values of /i, i.e., fi = 0.1 (left solid) and fi — 0.13 (right dashed); 
the corresponding labels are a) and b) for the symmetric (in the linear limit) branch shown above, while c) and d) are for the 
antisymmetric (in the linear limit) branch shown below. 
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FIG. 4: The value of fi at which branches III and IV disappear as a function of e with respect to Fig. (3] The blue solid 
lines and the green dashed-dotted lines denote the numerically found solutions and the result of the two-mode approximation, 
respectively. 
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to the symmetric ground eigenstate uq. The asymmetry arises immediately after the deviation from the Hnear hmit 
of — *■ 0, and becomes increasing as one drifts further away. Similarly, branch II starts from fi — toi, the eigenvalue 
of the anti-symmetric first excited eigenstate, and becomes increasingly asymmetric as it gets away from its linear 
limit, i.e., when N gets larger. The four bottom panels of Fig. [3] provide specific examples of the profiles of this 
continuation. 

Branches III and IV correspond to a pair of two other asymmetric solutions which collide at some critical value 
oi fi = fi'^'^ and disappear from then on through a saddle-mode bifurcation. As for the two examples in Fig. [31 the 
critical points are 0.1140 and 0.1054, for the cases e = 0.25 and e — 0.5, respectively. Moreover, branch IV is observed 
to move towards branch I as e — s- 0. Specifically, when e = 0, branch IV merges into branch I and the diagram turns 
into a pitchfork bifurcation (the two asymmetric branches are mirror images of each other and both bifurcate from the 
symmetric solution in that case), which is the case with attractive interactions analyzed in Ref. [2^ . As e increases, 
the critical point /i^'' of the saddle-node bifurcation decreases and tends to negative infinity rapidly, especially as e 
gets closer to 1; then, branches III and IV keep moving to the left of the diagram and disappear finally when e is close 
enough to 1. Also, the stability analysis indicates that branch III is the only unstable one among the four solutions, 
while the other three are all stable for any value of /i. 

As a stringent test of our analytical two-mode approximation, in Fig. |3]we show the critical point of the saddle-node 
bifurcation /^^'' as a function of the coUisional inhomogeneity parameter e. We observe that the solid line of the fully 
numerical results is almost identical to the dashed line yielding the theoretical prediction for the occurrence of this 
bifurcation. This figure demonstrates excellent agreement between the analytical results and the numerical findings. 




0.1 0.15 0.2 

M. 

FIG. 5: The norm of the solutions of Eq. ([4]) for a nonlinearity coefficient r(a;) with parameters such that a + /3 = — 1 and 
a — P — —1 + e, with e — 1.05 (top left), e = 1.5 (top right) and e — 2 (bottom), as a function of /i. The notation is the same 
as in Fig. [S] 

Next, we study the behavior of the solutions when e > 1, in which case the nonlinearity is no longer purely 
attractive. Remarkably, in this case, only two solutions, branch I and II, still survive for all the examined values of 
£ > 1. We realize that when e grows larger from to 1, both branches I and II move 'clockwise' and branch II appears 
to be almost vertical when e passes the value 1. As £ continues to increase from 1 to 2, we find an interesting new 
phenomenon: within a certain small range of iV, the solution II only exists for chemical potentials fi slightly less than 
the eigenvalue uji (corresponding to the eigenstate ui) before it meets a turning point; the latter occurs, for example, 
at /Lt = 0.1520 or 0.1554 for the cases e = 1.05 or 1.5, respectively. After the turning point, the solution exists when 
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/i gets larger. This phenomenon is shown in Fig. [5l in which it is observed that branch II starts at the linear limit, 
persists with dN/dfi < for a narrow interval of chemical potentials, before turning to the right with and acquiring 
dN/d^ > 0. It is also worth mentioning that both branches I and II are still asymmetric and stable in this case; thus, 
here we observe an interesting deviation from the well-known Vakhitov-Kolokolov criterion [59] (see, e.g., a relevant 
discussion in Ref. _.60]) about the slope of the branch determining its linear stability. We present an explanation 
of the relevant feature at the end of this section (after observing similar features in the principally defocusing case 
below). 

Now let us consider the case in which the parameters involved in the nonlinearity coefficient r(a;) are such that 

a + /3=l-e, a-/3 = l, (18) 
with < e < 2, as shown in Fig. [6l Notice that Eqs. (|12p and p8)) produce the same results when e = 2. 
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FIG. 6: Similar to Fig. [T] but with a nonlinearity coefficient F(a;) such that a + /3 — 1— e,a — P = l and b = 1, with e — 0.5 
(blue solid line), e = 1 (green dashed line) and e = 2 (red dashed-dotted line). 
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FIG. 7: The normalized number of particles A'^ for the solutions of Eq. (|4]) as a function of fi in the case of repulsive interatomic 
interactions, namely, for a nonlinearity coefficient F(a;) with parameters such that a + 13 — 1 — e and a — /3 = 1, with e = 0.25 
(left) and e = 0.5 (right). The notation is the same as in Fig. [3] 

Figure [7] shows some prototypical examples of the bifurcation diagram of the relevant solutions when < e < 1 , 
i.e., in the case of a purely repulsive nonlinearity (notice that the branches turn to the opposite direction than in Fig. 
[3|). Similar to the previous case, four solutions are found, two of which disappear through a saddle node bifurcation, 
allowing only two to survive in the non-interacting (linear) limit of N —^ 0. As before, the saddle-node bifurcation 
becomes a pitchfork one (but now emerging from the anti-symmetric branch) in the case of a coUisionally homogeneous 
environment i.e., for e — > 0. Preserving the original notations, branches I and II are the two asymmetric solutions 
extending to the linear limit, and starting at /x = Wq and ^ = wi, respectively. Once again branches III and IV 
disappear when decreases to some critical value ^'^^ , which goes up to infinity as e increases to 1. The dependence 
of fi'^^ on the inhomogeneity parameter e is shown in Fig. [51 It is important once again to highlight the good 
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qualitative and even quantitative (apart from the case of very large s) agreement of the two-mode prediction for /i' 
with the full numerical results. 
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FIG. 8: The critical value ^'^^ of the normalized chemical potential at which branches III and IV disappear as a function of e 
with respect to Fig. (6] The blue solid lines and the green dashed-dotted lines denote the numerically found solutions and the 
prediction of the two-mode approximation, respectively. 
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FIG. 9: The normalized number of particles of the solutions of Eq. Q as a function of fi, for a nonlinearity coefficient T{x) 
with parameters such that a + P = 1 — e and a — l3 — 1, with e — 1.05 (left) and e — 1.5 (right), The notation is the same as 
in Fig. H 



It is interesting to note that in this case as e increases and one of the wells becomes less repulsive (and eventually 
attractive for e > 1), the branches I and II keep rotating 'counterclockwise'. For e > 1, the turning of the branch arises 
again, on branch I this time, as shown in Fig. [9l In this case, the bifurcation diagram contains only two branches 
which, for sufficiently large e, feature opposite monotonicity of the dependence of N on the chemical potential /i 
(although both branches are linearly stable). Similarly to the previous case, branch III is the unstable solution, while 
the other three remain stable. In all the cases (and even these of large e), we again note the strong agreement between 
the bifurcation diagram predicted by the two-mode approximation, in comparison with the numerical results. 

We observe that in Fig. \5\ (e.g., in its top left panel) for the case where the nonlinearity is principally focusing, as 
well as in Fig. [n](e.g., in its left panel), where it is principally defocusing, one of the branches changes its monotonicity, 
as a result of the spatially dependent nonlinearity. As indicated previously, given the slope condition of the Vakhitov- 
Kolokolov criterion [59| , it appears to be rather surprising that this change of monotonicity is not accompanied by a 
change of stability. However, we argue here that it is not. Defining the well known linearization operators 

L+ = -^dl + V{x)+3g{x)\u\'-fi (19) 

L- = -^dl + V{x)+g{x)\u\' -f,, (20) 

it is known from the work of (6l| that when |n(L+) — n(L_)| — 1 [i.e., the number of negative eigenvalues of 
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minus the negative eigenvalues of L_ is in absolute value equal to 1], instability arises when the slope condition of 
Vakhitov-Kolokolov is violated. Since, when |n(L_|_) — n(L_)| > 1, the relevant theory indicates that the solution 
is always unstable, we argue that what should be happening here is that the violation of the slope condition of 
Vakhitov-Kolokolov is not associated with an instability because it occurs at the same time as the change of the 
count of — n(L_)| (from 1 to or vice versa). This is precisely what we numerically illustrate in Fig. [TOl In 

particular, the left panel concerns the principally focusing case associated to the first excited state "turning branch" 
of Fig. O There, given the fact that n(L^) = 1 (due to the one zero-crossing of the configuration itself, which is an 
eigenfunction of L_ with a eigenvalue), we expect the count of eigenvalues n(L+) to change from 2 to 1 exactly at 
the critical point, precisely as observed in the figure. On the other hand, the right panel concerning the ground state 
branch of Fig. [5] for the principally defocusing case, corresponds to a case with n{L^) = 0. Hence, when the slope 
condition is violated (near the linear limit in this case), it has to be true that n(L_|_) = 0, while after the turning 
point of the branch, the slope condition is satisfied and then n{L+) = 1, which again is consonant with the observed 
stability. These features are clearly illustrated in the right panel of Fig. [TOl We believe that similar considerations 
and counts may resolve apparent paradoxes that seem to be encountered e.g. in [60Ll6l|. 




FIG. 10: The eigenvalues of the operator L+ — + ^i^) + 3(?|mP — At, which involve negative parts, as a function of N. 

The left and right panels relate to the two "turning branches", i.e. branch II in Fig. [5] and branch I in Fig. |9] respectively, for 
the case that e = 1.05. 



IV. TWO-MODE DYNAMICS 



So far, we have considered the two-mode reduction as a tool for identifying (quite successfully, as shown above) the 
stationary states of the underlying problem. However, here we illustrate how the same tool can be used to understand 
the system dynamics. The two principal dynamical features of the (symmetric) double well system involve the 
oscillations of matter between the two wells (for low population asymmetries between the wells) and the nonlinearly 
induced self-trapping regime (for high population asymmetries between the wells); see e.g., [lEIill- For this reason, 
although our setting here is inherently asymmetric (due to the nature of r(x)), it is of interest to develop a variant 
of the two-mode approximation that accounts for the population imbalance between the wells. 

In order to formulate the problem based on the populations of the left and right well, one can reformulate our 
two-mode decomposition as 

u{x,t) = CR{t)^{uQ{x) + ui{x)) + CL{t)^{uo{x) - ui{x)), (21) 

where it is clear that cl and cn are connected with cq and ci through a simple linear transformation. If we then 
decompose ct — pie^'^' i — R, L and define the phase difference Acj) — (pL — <Pr, and the population imbalance 
z = pj^ — p\ (recall that N = p\ + p^), one can then project the equation to [uq + ui), as well as to (ug — ui) and 
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eventually obtain the dynamical equations for the conjugate variables z and A(/). These read: 

A,^ = --(ylo-lOS + Ai)+iV(i?o + ^i)-cos(A0)^=_=[/io-Mi + Y(^o-Ai)j 

_ z 
+ cos(A0)^==(i?o-i?i) -cos(2A(/.)-(Ao-2S + Ai) (22) 
V A*-^ — 4 

/ / N \ N"^ — z^ 

z= V^V23^sin(A?i)(^/xo-Aii + y(Ao-Ai)+z(Bo-Si)j + sin(2A0)(Ao - 2B + Ai) (23) 

It is interesting to examine the dynamical evolution of these equations (their stationary states are identical to the ones 
identified above), and to compare their dynamics with the corresponding PDE dynamics. In the latter, one can also 
make similar projections to (mq + ui), as well as to (mq — ui) and define accordingly cl.Bi ^ well as thereafter z and 
A0. A comparison of the (z, A(/i) phase space for different e of such canonically conjugate variables can be found in 
Fig. 111! It can be clearly observed that for non-zero e the reflection symmetry around z = is broken. For increasing 
e the population imbalance of the stationary states shifts towards larger values of z (see also Fig. [T^ . denoting a 
larger occupation of the right well compared to the left one. This can be understood by taking into account that 
we consider the case of repulsive interaction. There, an increase of e implies a reduction of interaction in the right 
well (see Fig. So, pictorially speaking, the atoms feel less repulsion in the right well than in the left one, leading 
to the observed population imbalance. For small populational asymmetries (around the nonzero stationary ones), 
the system executes inter-well matter-transfer oscillations. On the other hand, beyond a critical asymmetry in the 
initial populations, we enter the self-trapping regime, as is shown both in the partial differential equation of the full 
dynamics and in the above reduced two-mode description of Eqs. For non-zero e the value of this critical 

points taken at (50 = depend also on their sign, in contrary to the symmetric case for e = 0. For small values of e 
there are two critical points, which increase with increasing e. Eventually, the positive one ceases to exist as one can 
see by inspecting Fig. [11] for e = 0.5. 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 

A 0[;u] A (])[7t] 

FIG. 11: Phase diagrams for e = 0.25 (left column) and e = 0.5 (right column) and A'' = 0.5. The results from the ordinary 
differential equations (|22|l - (|23p are shown in the top row, while the corresponding results from the partial differential equation 
are shown in the bottom row. 

A detailed discussion of the dynamical comparison of the two-mode approximation with the corresponding results 
of the full system can be found in [l^. There, it is highlighted how the effect of a reduced barrier height of the 
trapping potential or increased interactions can deteriorate the dynamical effectiveness of the two-mode reduction. 

Here we explore the role of the collisional inhomogeneity in affecting the accuracy of this reduction, by illustrating a 
prototypical diagnostic, namely the critical population imbalance threshold beyond which self-trapping occurs, as this 
is identified in the full dynamical equation and as it is obtained within the two-mode reduction. This is shown in Fig. 
[T2] It can be seen that as e increases, both critical points (the positive and the negative one) increase. The positive 
one vanishes when it is equal to the total number of atoms. By comparing the results obtained by solving the PDE 
with the results of the ODEs one observes that the two-mode approximation becomes less accurate, as e increases, in 
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TDGPE 
fixpoints 
2M 



0.1 0.2 0.3 0.4 0.5 

e 

FIG. 12: The critical points Zcr of population imbalance above which there exists self-trapping are shown in the time-dependent 
Gross-Pitaevskii equation (TDGPE) and in the two-mode approximation (2M) for A'' = 0.5. Notice that in this setting, the 
population imbalance becomes nonzero when e 7^ due to the asymmetry (as shown by the dashed line). 



capturing the corresponding threshold, a feature that we have also seen regarding stationary state properties; see e.g. 
Fig. l 



V. CONCLUSIONS AND OUTLOOK 

In this paper, we studied in detail the nature of the most fundamental matter-waves (emanating from the linear 
limit of the problem) that emerge from a quasi- ID coUisionally inhomogeneous Bose-Einstein condensate trapped in 
a double-well potential. We specifically considered a setup which features distinct scattering lengths between the 
two wells, including also the case where the nonlinearity coefficient of the pertinent Gross-Pitaevskii equation has 
different signs in the the two wells. We observed that the relevant phenomenology is different in the considered 
coUisionally inhomogeneous environment in comparison to the coUisionally homogeneous case studied previously. In 
particular, even for weak inhomogeneities, the asymmetry (that the spatial dependence of the nonlinearity introduces) 
induces a modification of the bifurcation picture of the double-well potential and a change in the nature of the 
symmetry-breaking bifurcation from a pitchfork to a saddle-node; this is reminiscent of similar modifications to the 
bifurcation picture due to asymmetries of the linear potential [22]. On the other hand, it was found that as the 
strength of the inhomogeneity is increased, even this saddle-node "recollection" of the symmetry-breaking bifurcation 
eventually disappears (the corresponding critical point is pushed to infinity) and only one nonlinear branch persists 
that corresponds to each state of the problem's linear limit. Interestingly also, as the inhomogeneity acquires opposite 
sign values of the scattering length, the monotonicity of the number of atoms' dependence on the chemical potential 
may change (and may even be different between the two branches), although they do maintain their stability (which 
presents an interesting deviation -in this spatially inhomogeneous case- from the well-known Vakhitov-Kolokolov 
criterion that we rationalized in detail above). All of this phenomenology, including the detailed bifurcation diagram 
and other specific features, such as the critical point of the saddle-node bifurcation and its dependence on the degree 
of inhomogeneity, are captured remarkably accurately by a Galerkin-type, two-mode approximation and a resulting 
simple set of algebraic equations. Finally, we also briefly discussed dynamical aspects of the system including a phase 
space reduction in the population imbalance-interwell phase difference space. We have illustrated that despite the 
ensuing asymmetry in the phase space (for positive vs. negative population imbalances), the two- mode reduction 
can accurately capture the threshold of transition from Josephson matter-wave oscillations to nonlinearly induced 
self-trapping (especially so for weak collisional inhomogeneities). 

Our investigation herein presents some testable predictions for the original physical system and its realization in 
atomic physics or, also possibly, in nonlinear optics |25j ; in particular, as the nonlinearity becomes increasingly different 
in the two wells, the emergence of additional states should be occurring for increasingly larger number of atoms and 
eventually, such additional states (in particular, the stable one among them) should no longer be observable. It would 
be of particular interest if such spatial variations of the nonlinearity could be experimentally realized, which would 
allow the direct testing of our theoretical predictions. 
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